library(tidyverse)
library(broom)ETC1010/5510 Tutorial 11 - Solution
Introduction to Data Analysis
🎯 Workshop Objectives
- Introduction to linear models
- Interpreting linear models
- Linear model diagnostic
🔧 Instructions
Install the necessary pacakges
🌶️ Exercise: Introduction to Linear Model
Section A: Practice with Advertising Data
1. Read the data into RStudio. Remove the first column.
advert <- read_csv("data/Advertising.csv") %>%
select(-1)
head(advert)# A tibble: 6 × 4
TV radio newspaper sales
<dbl> <dbl> <dbl> <dbl>
1 230. 37.8 69.2 22.1
2 44.5 39.3 45.1 10.4
3 17.2 45.9 69.3 9.3
4 152. 41.3 58.5 18.5
5 181. 10.8 58.4 12.9
6 8.7 48.9 75 7.2Before fitting a model, let’s do some EDA. Plot the data with sales on the y-axis, and the rest of the variables on the x-axis in 3 separate plots. Then combine them into one plot.
p1 <- ggplot(advert,
aes(x = radio,
y = sales)) +
geom_point(alpha = 0.2)
p2 <- ggplot(advert,
aes(x = newspaper,
y = sales)) +
geom_point(alpha = 0.2)
p3 <- ggplot(advert,
aes(x = TV,
y = sales)) +
geom_point(alpha = 0.2)
gridExtra::grid.arrange(p1, p2, p3, nrow = 1)
2. Fit a model with just radio as the x variable.
The simple linear regression is of the form:
\[y_i=\beta_0+\beta_1x_{i}+\varepsilon_i, ~~~ i=1, \dots, n\]
Suppose we are interested in the amount of variation in sales that is explained by the variation in radio.
advert_m1 <- lm(sales ~ radio, data = advert)
advert_m1
Call:
lm(formula = sales ~ radio, data = advert)
Coefficients:
(Intercept) radio
9.3116 0.2025 The formula response ~ explanatory specifies the response variable and explanatory variable from the data. In this case, the response is sales in in thousands of units and the explanatory variable is spending on radio advertisement, in thousands of dollars.
Note that by default, the lm() function will include the intercept term.
We can extract components of the regression line using the broom package. The function tidy() will extract the estimated coefficients as a tidy table.
library(broom)
coef_advert_radio <- tidy(advert_m1)
coef_advert_radio# A tibble: 2 × 5
term estimate std.error statistic p.value
<chr> <dbl> <dbl> <dbl> <dbl>
1 (Intercept) 9.31 0.563 16.5 3.56e-39
2 radio 0.202 0.0204 9.92 4.35e-19The fitted model has the following functional form:
\[\widehat{\text{sales}} = 9.311 + 0.202~\text{radio}\]
The coefficients can be interpreted as
Slope: For each additional thousand dollars spending on radio advertising, on average, the sales increases by 202 units.
Intercept: When there is no advertising spending on radio, sales increases by 9 thousand units.
We can draw our fitted linear model with geom_smooth(method = "lm")
ggplot(advert,
aes(x = radio,
y = sales)) +
geom_point(alpha = 0.2) +
geom_smooth(method = "lm")
3. Diagnostics and augment
The residuals from the model are the differences between the observed values and the fitted values. The residuals and fitted values can be extracted from using augment() function as part of the broom package:
advert_m1_diagnostics <- augment(advert_m1)
advert_m1_diagnostics# A tibble: 200 × 8
sales radio .fitted .resid .hat .sigma .cooksd .std.resid
<dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 22.1 37.8 17.0 5.13 0.00982 4.27 0.00722 1.21
2 10.4 39.3 17.3 -6.87 0.0109 4.26 0.0143 -1.62
3 9.3 45.9 18.6 -9.31 0.0167 4.23 0.0409 -2.20
4 18.5 41.3 17.7 0.825 0.0124 4.29 0.000237 0.194
5 12.9 10.8 11.5 1.40 0.00854 4.28 0.000467 0.329
6 7.2 48.9 19.2 -12.0 0.0200 4.20 0.0822 -2.84
7 11.8 32.8 16.0 -4.15 0.00707 4.28 0.00339 -0.975
8 13.2 19.6 13.3 -0.0806 0.00531 4.29 0.000000952 -0.0189
9 4.8 2.1 9.74 -4.94 0.0152 4.27 0.0105 -1.16
10 10.6 2.6 9.84 0.762 0.0147 4.29 0.000241 0.180
# ℹ 190 more rowsWe can visualise the residuals directly on the scatter plot:
ggplot(advert,
aes(x = radio,
y = sales)) +
geom_point(alpha = 0.2)+
geom_smooth(method = "lm", se = FALSE) +
# overlay fitted values
geom_point(data = advert_m1_diagnostics,
aes(y = .fitted),
color = "blue",
alpha = 0.2) +
# draw a line segment from the fitted value to observed value
geom_segment(data = advert_m1_diagnostics,
aes(xend = radio, y = .fitted, yend = sales),
color = "blue",
alpha = 0.2)
Another way of visualising residuals is to produce a scatterplot of the residuals against the fitted values:
ggplot(advert_m1_diagnostics,
aes(x = .fitted, y = .resid)) +
geom_point()
Yet, another and a better way to visualise residuals is using resid_panel() which gives you a panel of diagnostic residual plots.
library(ggResidpanel)
resid_panel(advert_m1, plots = "all")
Recall the assumptions of the linear model from the lecture. What are they?
- Are they met from the plots?
- Are there observations with large positive residuals? What do they mean?
4. How good is your fit?
- \(R^2\) is a common measurement of strength of linear model fit.
- \(R^2\) tells us % variability in response explained by the model.
We can extract model fit summaries using glance()
advert_m1_summary <- glance(advert_m1)
advert_m1_summary# A tibble: 1 × 12
r.squared adj.r.squared sigma statistic p.value df logLik AIC BIC
<dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 0.332 0.329 4.27 98.4 4.35e-19 1 -573. 1153. 1163.
# ℹ 3 more variables: deviance <dbl>, df.residual <int>, nobs <int>The \(R^2\) value here can be interpreted as the proportion of variability in sales explained by spending on radio advertisement. That is approximately 33% of variability in sales is explained by spending on radio advertisement. Hmm, that is not much! 😒
Section B: Multiple Linear Regression
In reality, we will not only spend all the advertising spending on just one medium. We may choose to have other media of advertisement too such as newspaper and TV.
1. Run a multiple regression model to show the effect of different advertising media on sales. Obtain a summary of the regression results.
advert_m2 <- lm(sales~ TV + newspaper + radio,
data = advert)
summary(advert_m2)
Call:
lm(formula = sales ~ TV + newspaper + radio, data = advert)
Residuals:
Min 1Q Median 3Q Max
-8.8277 -0.8908 0.2418 1.1893 2.8292
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 2.938889 0.311908 9.422 <2e-16 ***
TV 0.045765 0.001395 32.809 <2e-16 ***
newspaper -0.001037 0.005871 -0.177 0.86
radio 0.188530 0.008611 21.893 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1.686 on 196 degrees of freedom
Multiple R-squared: 0.8972, Adjusted R-squared: 0.8956
F-statistic: 570.3 on 3 and 196 DF, p-value: < 2.2e-162. Is the relationship linear?
advert_m2_aug <- augment(advert_m2)
ggplot(data = advert_m2_aug,
aes(x = .fitted,
y = .resid)) +
geom_point() +
geom_hline(yintercept = 0, colour = "red")
3.Residual Plots
Residual plots can be used in order to identify non-linearity or other problems with our regression model. * Plot the residuals using resid_panel(). Take this out of the exercise, but leave in the solutions
library(ggResidpanel)
resid_panel(advert_m2, plots = "all")
Do you see any patterns? [ANS: Looks fine but the residuals tends to be more positive when the fitted values are at the lower or higher values.]
We can see that the residuals are not normal, in which the Q-Q plot, histogram and boxplot show data skewness.
The Cook’s Distance plot shows our data is having two extreme outliers. You may want to investigate on this.
👴 Exercise 9B: Data modelling using the gapminder data
Gapminder
- Hans Rosling was a Swedish doctor, academic and statistician, Professor of International Health at Karolinska Institute. Sadly he passed away in 2017.
- He developed a keen interest in health and wealth across the globe, and the relationship with other factors like agriculture, education, and energy.
- You can play with the gapminder data using animations at https://www.gapminder.org/tools/.
R package gapminder contains subset of the data on five year intervals from 1952 to 2007.
library(gapminder)
glimpse(gapminder)Rows: 1,704
Columns: 6
$ country <fct> "Afghanistan", "Afghanistan", "Afghanistan", "Afghanistan", …
$ continent <fct> Asia, Asia, Asia, Asia, Asia, Asia, Asia, Asia, Asia, Asia, …
$ year <int> 1952, 1957, 1962, 1967, 1972, 1977, 1982, 1987, 1992, 1997, …
$ lifeExp <dbl> 28.801, 30.332, 31.997, 34.020, 36.088, 38.438, 39.854, 40.8…
$ pop <int> 8425333, 9240934, 10267083, 11537966, 13079460, 14880372, 12…
$ gdpPercap <dbl> 779.4453, 820.8530, 853.1007, 836.1971, 739.9811, 786.1134, …Section A: Australia
1. How has life expectancy changed over the years, for each country?
Use a line plot to answer this question and provide your explanation here.
- There generally appears to be an increase in life expectancy in the different countries.
- A number of countries have big dips from the 1970s through 1990s.
- A cluster of countries starts off with low life expectancy but ends up closer to the highest life expectancy by the end of the period.
2. How has life expectancy in Australia changed over the years.
Use a line plot to answer this question and provide your explanation here.
oz <- gapminder %>%
filter(country == "Australia")
oz# A tibble: 12 × 6
country continent year lifeExp pop gdpPercap
<fct> <fct> <int> <dbl> <int> <dbl>
1 Australia Oceania 1952 69.1 8691212 10040.
2 Australia Oceania 1957 70.3 9712569 10950.
3 Australia Oceania 1962 70.9 10794968 12217.
4 Australia Oceania 1967 71.1 11872264 14526.
5 Australia Oceania 1972 71.9 13177000 16789.
6 Australia Oceania 1977 73.5 14074100 18334.
7 Australia Oceania 1982 74.7 15184200 19477.
8 Australia Oceania 1987 76.3 16257249 21889.
9 Australia Oceania 1992 77.6 17481977 23425.
10 Australia Oceania 1997 78.8 18565243 26998.
11 Australia Oceania 2002 80.4 19546792 30688.
12 Australia Oceania 2007 81.2 20434176 34435.3. Fit the linear model just for Australia: lifeExp ~ year
ggplot(data = oz,
aes(x = year,
y = lifeExp)) +
geom_line() 
4. Fit a linear model: lifeExp ~ year.
oz_lm <- lm(lifeExp ~ year, data = oz)
oz_lm
Call:
lm(formula = lifeExp ~ year, data = oz)
Coefficients:
(Intercept) year
-376.1163 0.2277 summary(oz_lm)
Call:
lm(formula = lifeExp ~ year, data = oz)
Residuals:
Min 1Q Median 3Q Max
-1.0250 -0.5201 0.1162 0.3781 0.7909
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -376.11630 20.54716 -18.30 5.09e-09 ***
year 0.22772 0.01038 21.94 8.67e-10 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.6206 on 10 degrees of freedom
Multiple R-squared: 0.9796, Adjusted R-squared: 0.9776
F-statistic: 481.3 on 1 and 10 DF, p-value: 8.667e-105. Tidy the model.
tidy(oz_lm)# A tibble: 2 × 5
term estimate std.error statistic p.value
<chr> <dbl> <dbl> <dbl> <dbl>
1 (Intercept) -376. 20.5 -18.3 5.09e- 9
2 year 0.228 0.0104 21.9 8.67e-10\[\widehat{lifeExp} = -376.1163 + 0.2277~year\]
6. Center year
Let us treat 1950 as the first year, so for model fitting, we are going to shift the years to begin in 1950. This makes interpretation easier.
gap <- gapminder %>%
mutate(year1950 = year - 1950)
oz <- gap %>%
filter(country == "Australia")7. Model for centered year
oz_lm <- lm(lifeExp ~ year1950, data = oz)
oz_lm
Call:
lm(formula = lifeExp ~ year1950, data = oz)
Coefficients:
(Intercept) year1950
67.9451 0.2277 8. Tidy the model centered year
tidy(oz_lm)# A tibble: 2 × 5
term estimate std.error statistic p.value
<chr> <dbl> <dbl> <dbl> <dbl>
1 (Intercept) 67.9 0.355 192. 3.70e-19
2 year1950 0.228 0.0104 21.9 8.67e-10\[\widehat{lifeExp} = 67.9 + 0.2277~year1950 \]
9. Extract residuals and fitted values using augment()
oz_aug <- augment(oz_lm, oz)
oz_aug# A tibble: 12 × 13
country continent year lifeExp pop gdpPercap year1950 .fitted .resid
<fct> <fct> <int> <dbl> <int> <dbl> <dbl> <dbl> <dbl>
1 Australia Oceania 1952 69.1 8691212 10040. 2 68.4 0.719
2 Australia Oceania 1957 70.3 9712569 10950. 7 69.5 0.791
3 Australia Oceania 1962 70.9 10794968 12217. 12 70.7 0.252
4 Australia Oceania 1967 71.1 11872264 14526. 17 71.8 -0.716
5 Australia Oceania 1972 71.9 13177000 16789. 22 73.0 -1.02
6 Australia Oceania 1977 73.5 14074100 18334. 27 74.1 -0.604
7 Australia Oceania 1982 74.7 15184200 19477. 32 75.2 -0.492
8 Australia Oceania 1987 76.3 16257249 21889. 37 76.4 -0.0508
9 Australia Oceania 1992 77.6 17481977 23425. 42 77.5 0.0505
10 Australia Oceania 1997 78.8 18565243 26998. 47 78.6 0.182
11 Australia Oceania 2002 80.4 19546792 30688. 52 79.8 0.583
12 Australia Oceania 2007 81.2 20434176 34435. 57 80.9 0.310
# ℹ 4 more variables: .hat <dbl>, .sigma <dbl>, .cooksd <dbl>, .std.resid <dbl>oz_aug$.fitted [1] 68.40051 69.53913 70.67775 71.81637 72.95499 74.09361 75.23223 76.37084
[9] 77.50946 78.64808 79.78670 80.92532oz_aug$.resid [1] 0.71948718 0.79086830 0.25224942 -0.71636946 -1.02498834 -0.60360723
[7] -0.49222611 -0.05084499 0.05053613 0.18191725 0.58329837 0.3096794910. Plot the year and lifeExp with points and add in the line.
ggplot(data = oz_aug,
aes(x = year,
y = .fitted)) +
geom_line(colour = "blue") +
geom_point(aes(x = year,
y = lifeExp))
11. Plot residuals against fitted values to reveal problems with the fit and interpret it.
ggplot(oz_aug,
aes(x = .fitted, y = .resid)) +
ylim(c(-5,5)) +
geom_point()
# Another way to look at Residuals: .std.resid with year:
ggplot(data = oz_aug,
aes(x = year,
y = .std.resid)) +
ylim(c(-5,5)) +
geom_hline(yintercept = 0,
colour = "white",
size = 2) +
geom_line() 
- Life expectancy has increased 2.3 years every decade, on average.
- There was a slow period from 1960 through to 1972, probably related to mortality during the Vietnam war.
Note: We use standardized residuals here is to make the residuals unitless. Therefore, it is easier for us to spot any unusual observations, which are those above the line of 2 and below -2. Remember the rules of 95% confidence interval?
Section B: Other countries??
1. Can we fit a model for New Zealand?
nz <- gap %>%
filter(country == "New Zealand")
nz_lm <- lm(lifeExp ~ year1950, data = nz)
nz_lm
Call:
lm(formula = lifeExp ~ year1950, data = nz)
Coefficients:
(Intercept) year1950
68.3013 0.1928 2. Can we fit a model for Japan?
japan <- gap %>%
filter(country == "Japan")
japan_lm <- lm(lifeExp ~ year1950, data = japan)
japan_lm
Call:
lm(formula = lifeExp ~ year1950, data = japan)
Coefficients:
(Intercept) year1950
64.4162 0.3529 3. Can we fit a model for Italy?
italy <- gap %>%
filter(country == "Italy")
italy_lm <- lm(lifeExp ~ year1950, data = italy)
italy_lm
Call:
lm(formula = lifeExp ~ year1950, data = italy)
Coefficients:
(Intercept) year1950
66.0574 0.2697 4. Is there a better way? What if we wanted to fit a model for ALL countries?
YES, let’s begin. 🙂
- Nest country level data (one row = one country)
by_country <- gap %>%
select(country, year1950, lifeExp, continent) %>%
group_by(country, continent) %>%
nest()
by_country# A tibble: 142 × 3
# Groups: country, continent [142]
country continent data
<fct> <fct> <list>
1 Afghanistan Asia <tibble [12 × 2]>
2 Albania Europe <tibble [12 × 2]>
3 Algeria Africa <tibble [12 × 2]>
4 Angola Africa <tibble [12 × 2]>
5 Argentina Americas <tibble [12 × 2]>
6 Australia Oceania <tibble [12 × 2]>
7 Austria Europe <tibble [12 × 2]>
8 Bahrain Asia <tibble [12 × 2]>
9 Bangladesh Asia <tibble [12 × 2]>
10 Belgium Europe <tibble [12 × 2]>
# ℹ 132 more rows- What is in the
datacolumn?
by_country$data[[1]]# A tibble: 12 × 2
year1950 lifeExp
<dbl> <dbl>
1 2 28.8
2 7 30.3
3 12 32.0
4 17 34.0
5 22 36.1
6 27 38.4
7 32 39.9
8 37 40.8
9 42 41.7
10 47 41.8
11 52 42.1
12 57 43.8It’s a tibble! 🤪
Let’s fit a linear model to each one of them.
lm_afganistan <- lm(lifeExp ~ year1950,
data = by_country$data[[1]])
lm_albania <- lm(lifeExp ~ year1950,
data = by_country$data[[2]])
lm_algeria <- lm(lifeExp ~ year1950,
data = by_country$data[[3]])🙈 But, we are copying and pasting this code more than twice… Is there a better way?
A case for map???
map(<data object>, <function>)
1. Write a function to fit lm to all countries and map the function to data column of by_country
fit_lm <- function(x){
lm(lifeExp ~ year1950, data = x)
}
mapped_lm <- map(by_country$data, fit_lm)
head(mapped_lm)[[1]]
Call:
lm(formula = lifeExp ~ year1950, data = x)
Coefficients:
(Intercept) year1950
29.3566 0.2753
[[2]]
Call:
lm(formula = lifeExp ~ year1950, data = x)
Coefficients:
(Intercept) year1950
58.5598 0.3347
[[3]]
Call:
lm(formula = lifeExp ~ year1950, data = x)
Coefficients:
(Intercept) year1950
42.2364 0.5693
[[4]]
Call:
lm(formula = lifeExp ~ year1950, data = x)
Coefficients:
(Intercept) year1950
31.7080 0.2093
[[5]]
Call:
lm(formula = lifeExp ~ year1950, data = x)
Coefficients:
(Intercept) year1950
62.2250 0.2317
[[6]]
Call:
lm(formula = lifeExp ~ year1950, data = x)
Coefficients:
(Intercept) year1950
67.9451 0.2277 2. Map inside the data?
country_model <- by_country %>%
mutate(model = map(data, function(x){
lm(lifeExp ~ year1950, data = x)
})
)
country_model# A tibble: 142 × 4
# Groups: country, continent [142]
country continent data model
<fct> <fct> <list> <list>
1 Afghanistan Asia <tibble [12 × 2]> <lm>
2 Albania Europe <tibble [12 × 2]> <lm>
3 Algeria Africa <tibble [12 × 2]> <lm>
4 Angola Africa <tibble [12 × 2]> <lm>
5 Argentina Americas <tibble [12 × 2]> <lm>
6 Australia Oceania <tibble [12 × 2]> <lm>
7 Austria Europe <tibble [12 × 2]> <lm>
8 Bahrain Asia <tibble [12 × 2]> <lm>
9 Bangladesh Asia <tibble [12 × 2]> <lm>
10 Belgium Europe <tibble [12 × 2]> <lm>
# ℹ 132 more rows3. Case for map (shorthand function)
country_model <- by_country %>%
mutate(model = map(data, ~lm(lifeExp ~ year1950,
data = .)))
country_model# A tibble: 142 × 4
# Groups: country, continent [142]
country continent data model
<fct> <fct> <list> <list>
1 Afghanistan Asia <tibble [12 × 2]> <lm>
2 Albania Europe <tibble [12 × 2]> <lm>
3 Algeria Africa <tibble [12 × 2]> <lm>
4 Angola Africa <tibble [12 × 2]> <lm>
5 Argentina Americas <tibble [12 × 2]> <lm>
6 Australia Oceania <tibble [12 × 2]> <lm>
7 Austria Europe <tibble [12 × 2]> <lm>
8 Bahrain Asia <tibble [12 × 2]> <lm>
9 Bangladesh Asia <tibble [12 × 2]> <lm>
10 Belgium Europe <tibble [12 × 2]> <lm>
# ℹ 132 more rows4. What’s the model?
country_model$model[[1]]
Call:
lm(formula = lifeExp ~ year1950, data = .)
Coefficients:
(Intercept) year1950
29.3566 0.2753 5. How do we summarise the content in the model?
tidy(country_model$model[[1]])# A tibble: 2 × 5
term estimate std.error statistic p.value
<chr> <dbl> <dbl> <dbl> <dbl>
1 (Intercept) 29.4 0.699 42.0 1.40e-12
2 year1950 0.275 0.0205 13.5 9.84e- 86. So should we repeat it for each one?
tidy(country_model$model[[1]])# A tibble: 2 × 5
term estimate std.error statistic p.value
<chr> <dbl> <dbl> <dbl> <dbl>
1 (Intercept) 29.4 0.699 42.0 1.40e-12
2 year1950 0.275 0.0205 13.5 9.84e- 8tidy(country_model$model[[2]])# A tibble: 2 × 5
term estimate std.error statistic p.value
<chr> <dbl> <dbl> <dbl> <dbl>
1 (Intercept) 58.6 1.13 51.7 1.79e-13
2 year1950 0.335 0.0332 10.1 1.46e- 6tidy(country_model$model[[3]])# A tibble: 2 × 5
term estimate std.error statistic p.value
<chr> <dbl> <dbl> <dbl> <dbl>
1 (Intercept) 42.2 0.756 55.8 8.22e-14
2 year1950 0.569 0.0221 25.7 1.81e-10NO!! There’s a better way. 😙
country_model %>%
mutate(tidy = map(model, tidy))# A tibble: 142 × 5
# Groups: country, continent [142]
country continent data model tidy
<fct> <fct> <list> <list> <list>
1 Afghanistan Asia <tibble [12 × 2]> <lm> <tibble [2 × 5]>
2 Albania Europe <tibble [12 × 2]> <lm> <tibble [2 × 5]>
3 Algeria Africa <tibble [12 × 2]> <lm> <tibble [2 × 5]>
4 Angola Africa <tibble [12 × 2]> <lm> <tibble [2 × 5]>
5 Argentina Americas <tibble [12 × 2]> <lm> <tibble [2 × 5]>
6 Australia Oceania <tibble [12 × 2]> <lm> <tibble [2 × 5]>
7 Austria Europe <tibble [12 × 2]> <lm> <tibble [2 × 5]>
8 Bahrain Asia <tibble [12 × 2]> <lm> <tibble [2 × 5]>
9 Bangladesh Asia <tibble [12 × 2]> <lm> <tibble [2 × 5]>
10 Belgium Europe <tibble [12 × 2]> <lm> <tibble [2 × 5]>
# ℹ 132 more rowsData Wrangle and tidy:
country_coefs <- country_model %>%
mutate(tidy = map(model, tidy)) %>%
unnest(tidy) %>%
select(country, continent, term, estimate)
country_coefs# A tibble: 284 × 4
# Groups: country, continent [142]
country continent term estimate
<fct> <fct> <chr> <dbl>
1 Afghanistan Asia (Intercept) 29.4
2 Afghanistan Asia year1950 0.275
3 Albania Europe (Intercept) 58.6
4 Albania Europe year1950 0.335
5 Algeria Africa (Intercept) 42.2
6 Algeria Africa year1950 0.569
7 Angola Africa (Intercept) 31.7
8 Angola Africa year1950 0.209
9 Argentina Americas (Intercept) 62.2
10 Argentina Americas year1950 0.232
# ℹ 274 more rowstidy_country_coefs <- country_coefs %>%
pivot_wider(id_cols = c(country, continent),
names_from = term,
values_from = estimate) %>%
rename(intercept = `(Intercept)`)
tidy_country_coefs# A tibble: 142 × 4
# Groups: country, continent [142]
country continent intercept year1950
<fct> <fct> <dbl> <dbl>
1 Afghanistan Asia 29.4 0.275
2 Albania Europe 58.6 0.335
3 Algeria Africa 42.2 0.569
4 Angola Africa 31.7 0.209
5 Argentina Americas 62.2 0.232
6 Australia Oceania 67.9 0.228
7 Austria Europe 66.0 0.242
8 Bahrain Asia 51.8 0.468
9 Bangladesh Asia 35.1 0.498
10 Belgium Europe 67.5 0.209
# ℹ 132 more rows7. Check for Australia
tidy_country_coefs %>%
filter(country == "Australia")# A tibble: 1 × 4
# Groups: country, continent [1]
country continent intercept year1950
<fct> <fct> <dbl> <dbl>
1 Australia Oceania 67.9 0.228Section C: Extension Exercise
- Fit the models to all countries.
- Pick your favourite country (other than Australia), print the coefficients, and make a hand sketch of the the model fit.
1. Plot all the models
country_aug <- country_model %>%
mutate(augmented = map(model, augment)) %>%
unnest(augmented)
country_aug# A tibble: 1,704 × 12
# Groups: country, continent [142]
country continent data model lifeExp year1950 .fitted .resid .hat
<fct> <fct> <list> <list> <dbl> <dbl> <dbl> <dbl> <dbl>
1 Afghanistan Asia <tibble> <lm> 28.8 2 29.9 -1.11 0.295
2 Afghanistan Asia <tibble> <lm> 30.3 7 31.3 -0.952 0.225
3 Afghanistan Asia <tibble> <lm> 32.0 12 32.7 -0.664 0.169
4 Afghanistan Asia <tibble> <lm> 34.0 17 34.0 -0.0172 0.127
5 Afghanistan Asia <tibble> <lm> 36.1 22 35.4 0.674 0.0991
6 Afghanistan Asia <tibble> <lm> 38.4 27 36.8 1.65 0.0851
7 Afghanistan Asia <tibble> <lm> 39.9 32 38.2 1.69 0.0851
8 Afghanistan Asia <tibble> <lm> 40.8 37 39.5 1.28 0.0991
9 Afghanistan Asia <tibble> <lm> 41.7 42 40.9 0.754 0.127
10 Afghanistan Asia <tibble> <lm> 41.8 47 42.3 -0.534 0.169
# ℹ 1,694 more rows
# ℹ 3 more variables: .sigma <dbl>, .cooksd <dbl>, .std.resid <dbl>p1 <- gapminder %>%
ggplot(aes(year, lifeExp, group = country)) +
geom_line(alpha = 1/3) + ggtitle("Data")
p2 <- ggplot(country_aug) +
geom_line(aes(x = year1950 + 1950,
y = .fitted,
group = country),
alpha = 1/3) +
xlab("year") +
ggtitle("Fitted models")library(gridExtra)
grid.arrange(p1, p2, ncol=2)
2. Plot all the model coefficients
p <- ggplot(tidy_country_coefs,
aes(x = intercept,
y = year1950,
colour = continent,
label = country)) +
geom_point(alpha = 0.5,
size = 2) +
scale_color_brewer(palette = "Dark2")library(plotly)
ggplotly(p)Let’s summarise the information we learnt from the model coefficients.
- Generally, the relationship is negative: this means that if a country starts with a high intercept, it tends to have lower rate of increase.
- There is a difference across the continents: Countries in Europe and Oceania tend to start with a higher life expectancy and increased; countries in Asia and America tended to start lower but have high rates of improvement; Africa tends to start lower and have a huge range in rate of change.
- Three countries had negative growth in life expectancy: Rwand, Zimbabwe, Zambia
Model diagnostics by country
country_glance <- country_model %>%
mutate(glance = map(model, glance)) %>%
unnest(glance)
country_glance# A tibble: 142 × 16
# Groups: country, continent [142]
country continent data model r.squared adj.r.squared sigma statistic
<fct> <fct> <list> <list> <dbl> <dbl> <dbl> <dbl>
1 Afghanistan Asia <tibble> <lm> 0.948 0.942 1.22 181.
2 Albania Europe <tibble> <lm> 0.911 0.902 1.98 102.
3 Algeria Africa <tibble> <lm> 0.985 0.984 1.32 662.
4 Angola Africa <tibble> <lm> 0.888 0.877 1.41 79.1
5 Argentina Americas <tibble> <lm> 0.996 0.995 0.292 2246.
6 Australia Oceania <tibble> <lm> 0.980 0.978 0.621 481.
7 Austria Europe <tibble> <lm> 0.992 0.991 0.407 1261.
8 Bahrain Asia <tibble> <lm> 0.967 0.963 1.64 291.
9 Bangladesh Asia <tibble> <lm> 0.989 0.988 0.977 930.
10 Belgium Europe <tibble> <lm> 0.995 0.994 0.293 1822.
# ℹ 132 more rows
# ℹ 8 more variables: p.value <dbl>, df <dbl>, logLik <dbl>, AIC <dbl>,
# BIC <dbl>, deviance <dbl>, df.residual <int>, nobs <int>3. Plot the \(R^2\) values using a histogram.
ggplot(country_glance,
aes(x = r.squared)) +
geom_histogram()
4. Countries with worst fit
Examine the countries with the worst fit, countries with \(R^2<0.45\), by making scatterplots of the data, with the linear model overlaid.
badfit <- country_glance %>% filter(r.squared <= 0.45)
gap_bad <- gap %>% filter(country %in% badfit$country)
gg_bad_fit <-
ggplot(data = gap_bad,
aes(x = year,
y = lifeExp)) +
geom_point() +
facet_wrap(~country) +
scale_x_continuous(breaks = seq(1950,2000,10),
labels = c("1950", "60","70", "80","90","2000")) +
geom_smooth(method = "lm",
se = FALSE)gg_bad_fit
Each of these countries was moving on a nice trajectory of increasing life expectancy, but later suffered a big dip during the time period.